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Abstract 

The  analysis  of  thin  structural  components,  which  are  characteristic  of  a  broad  class  of  Micro 
Air  Vehicles,  is  presented  herein.  A  direct  solution  approach  in  co-simulation  with  fluid-dynamics 
solvers  is  used.  An  original  variational  formulation  is  developed  for  the  inverse  problem  of  recon¬ 
structing  full- field  structural  displacement  and  pressure  distribution  of  membrane  wings  subjected 
to  static  and  unsteady  loads  from  membrane  strain  distribution.  Moving  Least  Squares  are  used  to 
smooth  and  remap  surface  strain  measurements,  estimated  from  Digital  Image  Correlation  (DIC), 
as  needed  by  the  inverse  solution  meshing.  The  same  approach  is  used  to  map  the  structural  and 
fluid  interface  kinematics  and  loads  during  the  fluid-structure  co-simulation.  The  inverse  analysis 
is  verified  by  reconstructing  the  deformed  solution  obtained  with  a  corresponding  direct  formu¬ 
lation,  based  on  nonlinear  membrane  structural  analysis  implemented  in  a  free  general-purpose 
multibody  dynamics  solver  and  tightly  coupled  in  co-simulation  with  a  CFD  solver.  Both  the 
direct  and  the  inverse  analyses  are  validated  by  comparing  the  direct  predictions  and  the  recon¬ 
structed  deformations  with  experimental  data  for  prestressed  rectangular  membranes  subjected 
to  static  and  unsteady  loads.  The  load  distributions  reconstructed  using  the  inverse  analysis  are 
compared  with  the  corresponding  ones  obtained  using  the  direct  analysis.  The  inverse  analysis 
runs  on  standard  off-the-shelf  PCs  and  can  be  implemented  in  real-time,  providing  load  distribu¬ 
tion  estimates  at  a  rate  in  the  order  of  tens  of  datasets  per  second. 
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Chapter  1 


Formulation 


1.1  Introduction 

The  work  presented  within  this  paper  seeks  to  obtain  full-field  estimations  of  the  structural 
displacement  field  for  a  membrane  wing  for  Micro  Air  Vehicle  (MAV)  applications.  The  source 
of  these  estimations  is  the  elastic  wing  deformation,  experimentally  measured  in  a  low-speed 
wind  tunnel  using  a  full-field,  non-contact,  digital  image  correlation  (DIC)  technique,  originally 
developed  by  researchers  at  the  University  of  South  Carolina  [30] . 

MAVs  were  defined  by  the  Defense  Advanced  Research  Projects  Agency  as  aircraft  with 
wingspan  less  than  15  cm  and  a  maximum  speed  less  than  15  m/s.  These  aircraft  can  be  utilized  for 
a  variety  of  missions,  carrying  payload  such  as  surveying  and  sensing  equipment.  In  MAV  designs, 
much  like  natural  fliers,  compliant  membranes  are  used  to  passively  enhance  flight  characteristics. 
Membrane  wings  display  unique  aerodynamic  characteristics  due  to  their  aeroelastic  nature,  which 
can  provide  performance  improvements  over  their  rigid  airfoil  counterparts.  Extensive  research 
has  been  conducted  on  the  dynamics  of  flexible  wings,  but  limited  research  exists  on  practical, 
computationally  predictive  models  dealing  with  the  dynamic  behavior  surrounding  these  wings. 

Stanford  et  al.  [28]  analyzed  the  effect  of  a  flexible  membrane  on  a  fixed  wing  with  both  ex¬ 
periments  and  numerical  modeling.  They  modeled  the  membrane  as  inextensible,  using  a  linear 
stress-stiffening  model.  The  linear  stress-strain  assumption  held  well,  because  the  strain  accumu¬ 
lated  due  to  the  aerodynamic  load  was  small  in  comparison  with  the  prestrain  of  the  membrane. 
Utilizing  experimentally  measured  deformations  of  a  membrane  wing  in  a  wind  tunnel  they  suc¬ 
cessfully  predicted  the  aerodynamic  forces  over  the  wing  for  a  range  of  flight  speeds  and  angles 
of  attack:  an  inverse  technique  was  formulated  by  locating  a  pressure  field  that  minimizes  the 
least-squares  difference  between  experimental  displacements  and  calculated  values;  the  calculated 
pressure  field  was  then  used  to  qualitatively  describe  the  important  flow  characteristic  over  the 
membrane  wing,  including  flow  stagnation,  pressure  recovery,  flow  reattachment  and  wing  tip 
vortices. 

Because  of  unsteadiness,  flow  separation,  and  turbulence,  previous  studies  using  panel  methods 
and  simplified  laminar  solvers  have  failed  to  capture  the  exact  effect  of  membrane  wings  on 
performance.  Gopalakrishnan  and  Tafti  [15]  addressed  this  issue  by  analyzing  flapping  flight  for 
a  flexible  wing  (at  low  Reynolds  number)  using  an  unsteady  large-eddy  simulation  (LES)  flow 
solver  coupled  with  a  linear  elastic  membrane  wing  model.  The  focus  of  the  study  was  to  evaluate 
the  effect  of  aeroelastic  cambering  on  flapping  flight  performance  using  a  linear  elastic  membrane 
model  with  different  prestress  values.  The  wing  is  treated  as  an  elastic  membrane  with  in-plane 
prestresses.  The  prestresses  of  the  wing  are  tailored  to  induce  a  camber  in  the  range  of  0.1-0.25 
times  the  chord  length,  and  their  effect  is  analyzed  based  on  changes  in  flow  structure  and  on 
variation  of  thrust  and  lift.  They  showed  that  the  introduction  of  camber  increases  thrust  and  lift 
production  significantly,  although  the  transverse  displacement  was  up  to  25%  of  the  chord,  which 
means  a  membrane  strain  of  more  than  12%,  which  in  turn  may  cause  substantial  change  of  the 
membrane  stress,  making  the  problem  strongly  nonlinear.  In  practice,  the  linearized  structural 
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Chapter  1 .  Formulation 


1.2.  Methodology 


model,  which  assumes  constant  membrane  prestress,  is  inadequate  for  straining  of  such  magnitude. 

Inverse  FEM  analysis  has  been  proposed  also  to  identify  the  mechanical  properties  of  mem¬ 
branes,  e.g.  when  for  nontraditional  methods  are  needed  for  characterizing  the  behavior  of  ma¬ 
terials  when  the  associated  boundary  problems  is  complex  [18,  17].  These  methods  finds  several 
applications  for  example  in  the  identification  of  the  mechanical  properties  of  biological  tissues  [1]. 

The  present  work  discusses  a  variational  principle  that  provides  the  mathematical  framework 
from  which  a  robust  inverse  finite  element  method  (IFEM)  is  developed.  The  problem  at  hand 
is  the  reconstruction  of  the  three-dimensional  deformations  of  membrane  structures  based  upon 
the  experimentally  measured  (discrete)  surface  strains  (and  well-defined  boundary  restraints), 
experienced  by  a  flexible  wing  during  flight,  and  the  estimation  of  the  aerodynamic  pressure 
exerted  on  the  wing.  With  an  estimate  of  the  pressure  distribution,  aerodynamic  loads  can  be 
estimated,  as  attempted  for  example  in  Carpenter  and  Albertani  [  0].  The  actual  loads  that 
cause  the  deformations  are  unknown;  however,  their  influences  are  represented  in  the  measured 
strains.  This  “inverse”  technique  (as  opposed  to  a  conventional  “direct”  technique:  estimating 
the  displacement  field  from  a  measured  pressure  distribution)  represents  a  viable  alternative  to 
conventional  pressure  measurement  techniques  in  low  Reynolds  number  environments.  In  fact, 
the  thin  elastic  membrane  wing  skins  used  to  decrease  the  vehicle  weight  and  obtain  a  certain 
amount  of  passive  shape  adaptation  [26]  are  particularly  susceptible  to  intrusive  measurements. 

Currently,  the  numerical  validation  of  the  flow  field  created  by  a  MAV  wing  is  largely  limited  to 
(i)  a  comparison  of  numerical  aerodynamic  coefficients  with  those  garnered  through  wind  tunnel 
test  analysis  [29],  or  (ii)  a  comparison  with  flow  visualization,  focusing  on  the  flow  separation, 
transition,  and  reattachment  locations  over  the  wing  [20].  Knowledge  of  the  full-field  differential 
pressure  distribution  over  the  wing  surface  can  provide  a  further  level  of  comparison,  indicating 
areas  over  the  wing  where  the  model  may  be  inadequate.  An  inverse  method  could  take  the 
deformed  wing  shapes  and  estimate  the  resulting  pressure  distribution  within  flight  regimes  that 
are  difficult,  if  not  impossible  to  simulate  through  either  CFD  or  wind  tunnel  testing. 

Aerodynamically,  inverse  problems  have  two  main  applications:  they  could  be  used  for  inverse 
design  problems  for  optimal  airfoil  geometries  [14] ,  or  for  structural  health  monitoring  (an  elastic 
wing  is  mounted  with  deformation  sensors,  typically  strain  gages  or  fiber-optics,  whose  signals 
are  used  to  reconstruct  the  displacement  field  [31],  or  the  original  wing  loading  [25]). 

1.2  Methodology 

A  membrane  structural  model  is  developed  for  both  direct  and  inverse  dynamics.  The  direct 
dynamics  analysis  is  used  to  predict  the  deformed  shape  under  specified  loads.  The  inverse  kine¬ 
matics  analysis  is  used  to  reconstruct  the  membrane  shape  from  the  membrane  strain  field.  The 
inverse  dynamics  analysis  is  used  to  reconstruct  the  pressure  distribution. 

A  membrane  finite  element,  implemented  in  a  multibody  formulation  [22],  is  used  in  co¬ 
simulation  with  a  fluid  dynamics  solver  to  predict  the  configuration  of  the  system  under  static 
and  unsteady  loads  [2,  3]. 

An  original  approach  is  developed  for  the  inverse  problem  of  the  reconstruction  of  full-field 
structural  displacements  of  membrane  wings  utilizing  surface  strain  measurements  [4,  6]. 

The  inverse  problem  of  full-field  structural  displacement  reconstruction  is  addressed  through 
the  application  of  a  variational  formulation,  leading  to  a  versatile,  robust  and  computationally 
efficient  inverse  membrane  nonlinear  finite  element  analysis  [5],  which  was  inspired  by  analogous, 
although  linear,  approaches  developed  in  the  past  for  shell-like  structures,  see  Shkarayev  et  al. 
[25],  Tessler  and  Spangler  [31]. 

In  the  current  case,  nonlinear  elasticity  is  mandatory  to  capture  the  essence  of  the  transverse 
load  carrying  capability  of  membranes,  whereas  in  the  previous  mentioned  prior  formulations 
the  problem  was  restricted  to  linear  elasticity:  when  subjected  to  a  finite  amount  of  transverse 
displacement,  the  assumption  of  constant  membrane  prestress  used  in  linearized  membrane  models 
is  no  longer  acceptable.  The  complete  set  of  membrane  strain  measures,  consistent  with  non-linear 
membrane  theory,  need  to  be  used. 
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Chapter  1.  Formulation 


1.3.  Origin  of  Least-Squares  Problems 


Exploiting  the  functionalities  provided  by  the  free  software  project  FEniCS1  (a  collection  of 
libraries  specifically  designed  for  the  automated  and  efficient  solution  of  PDEs),  a  three- node 
inverse  membrane  element  was  developed  (see  Alioli  et  al.  [4]):  three  displacement  degrees  of 
freedom  are  used  for  each  node,  namely  two  displacement  components  in  the  plane  of  the  mem¬ 
brane  and  one  along  the  transverse  direction.  The  error  function  is  the  difference  between  the 
membrane  strain  measures  expressed  as  functions  of  the  displacements  and  the  corresponding 
membrane  strain  measures  obtained  from  the  experimental  strains  by  re-sampling. 

A  penalty-parameter  controlled  regularization  term  mitigates  the  ill-posedness  of  the  problem 
associated  with  the  non-uniqueness  of  the  solution  in  terms  of  transverse  displacement  for  given 
membrane  stresses  and  with  the  high-order  nonlinearity  of  the  membrane  strains  with  respect  to 
transverse  displacement.  In  fact,  in  addition  to  the  usual  level  of  ill-posedness  of  linear  inverse 
problems  (they  do  not  necessarily  satisfy  conditions  of  existence,  uniqueness,  and  stability,  see 
for  example  Bakushinsky  and  Goncharsky  [9],  Shkarayev  et  al.  [25]),  the  present  one  is  also 
characterized  by  the  fact  that  for  null  or  low  membrane  prestress  the  problem  is  exactly  singular 
in  configurations  that  present  no  transverse  displacement  of  the  membrane. 

The  reconstructed  shape  of  the  membrane  is  used  to  estimate  the  surface  loads.  The  procedure 
is  verified  and  validated  by  correlation  with  the  surface  load  values  predicted  by  the  coupled  fluid- 
structure  analysis. 

The  present  work  uses  an  experimental  setup  that  can  accurately  obtain  the  full-field  three- 
dimensional  displacement  and  membrane  strain  over  a  moderate  size  wing  in  wind  tunnel  testing 
conditions.  The  proposed  methodology  enables  accurate  reconstruction  of  the  three-dimensional 
displacement  field.  It  may  be  effectively  employed  to  develop  real-time  processing  of  the  sensed 
information. 

Analytical  and  numerical  results,  along  with  experimental  measurements  of  actual  membrane 
wing  artifacts  subjected  to  a  variety  of  steady  and  unsteady  flow  conditions  are  used  to  validate 
the  proposed  formulation. 

Experimental  data  is  based  on  DIC  in  conjunction  with  a  load  cell  and  tensile  test  frame  to 
measure  stress  and  strains:  DIC  measurements  were  taken  to  generate  virtual  strain  sensors  on  the 
surface  of  the  membrane  [10].  Measurements  are  further  manipulated  using  moving  least-squares 
(MLS)  [24]  to  remap  the  measured  displacements  and  strains  on  the  same  grid  that  is  used  for 
the  inverse  analysis. 

Historically,  due  to  its  commercial  availability,  the  membrane  used  in  the  experiments  has 
been  made  of  isotropic  material.  Methods  such  as  pretension  before  mounting  it  on  a  frame 
have  been  developed  to  alter  its  characteristics.  However,  by  developing  a  non-isotropic  elastic 
membrane  material  capable  of  being  tailored  to  an  applicable  stiffness  range  for  membrane- 
based  wings,  researchers  might  be  able  to  better  replicate  the  successful  characteristics  of  natural 
fliers.  A  non-isotropic  membrane,  in  fact,  should  response  differently  to  pressure  loading  between 
the  longitudinal  and  transverse  direction:  this  response  could  be  used  to  vary  designed  flight 
characteristics  of  future  membrane  wings  [32].  Thus,  a  hydrostatic  membrane  pressure  test  was 
also  conducted  to  characterize  the  behavior  of  the  non-isotropic  membrane  under  a  constant  and 
uniform  pressure  distribution.  A  non-isotropic  and  a  silicone  control  sample  were  each  secured  over 
a  frame  and  subjected  to  a  pressure  differential  across  the  membrane.  The  membrane  response 
to  the  pressure  differential  is  used  to  predict  the  membrane  behavior  under  aerodynamic  loading. 
As  researchers  refine  membrane  wing  designs,  the  availability  of  elastic,  non-isotropic  material 
similar  to  bat  wings  will  provide  another  tool  in  the  development  of  a  functional  MAV  membrane 
design  with  potential  for  expanding  their  flight  envelope. 


1.3  Origin  of  Least-Squares  Problems 

In  this  section,  the  basics  of  the  Least-Squares  problem  are  recalled,  to  provide  the  essential 
background  for  the  formulation  of  the  inverse  problem  that  is  developed  in  Section  1.4. 

xhttp : //f eniesproj  ect . org/ 


Real-Time  Pressure  Distribution  Estimation  on  Wings  Via  Displacements  and  Strains 


4 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


Chapter  1.  Formulation 


1.3.  Origin  of  Least-Squares  Problems 


A  variety  of  practical  problems  can  be  formulated  as  the  minimization  of  the  square  of  an  error 
function;  i.e.,  of  a  scalar  function  F(x)  which  can  be  expressed  as  a  sum  of  squares  of  nonlinear 
functions: 

F(x)=2^/,(x)2=2lf(x)f  (1.1) 

2=1 

where  x  =  (xi,x2,  ■  ■  ■  ,xn)  is  a  vector  and  each  /;  is  a  function  from  Rn  to  R.  The  factor  t  has  been 
included  in  (1.1)  in  order  to  avoid  the  appearance  of  a  factor  2  in  the  derivatives.  The  fi  are 
referred  to  as  residuals.  It  is  assumed  that  m>n. 

Although  function  (1.1)  can  be  minimized  by  a  variety  of  unconstrained  methods,  in  most 
circumstance  the  properties  of  (1.1)  make  it  worthwhile  to  use  methods  designed  specifically  for 
the  least-squares  problem.  In  particular,  the  gradient  and  Hessian  of  (1.1)  have  a  special  structure. 
Let  the  m  x  n  Jacobian  matrix  of  f(x)  be  denoted  by  J(x),  and  let  matrix  Gj(x)  denote  the  Hessian 
matrix  of  /;(x).  Then 

g(x)  =  J(x)Tf(x)  (1.2) 

G(x)  =  J(x)TJ(x)  +  Q(x)  (1-3) 

where  Q(x)  =  J27L i  /;(x)Gk(x)-  From  (1.3)  we  can  observe  that  the  Hessian  of  a  least-squares  objec¬ 
tive  function  consists  of  a  special  combination  of  first-  and  second-order  information.  Typically, 
least-squares  methods  are  based  on  the  hypothesis  that  eventually  the  first-order  term  J(x)TJ(x) 
of  (1.3)  becomes  dominant  with  respect  to  the  second-order  term  Q(x).  This  assumption  is  not 
justified  when  the  residuals  at  the  solution  are  very  large:  in  such  a  case,  one  might  as  well  use  a 
general  unconstrained  method.  However,  for  many  problems,  the  residual  at  the  solution  is  small 
enough  to  justify  the  use  of  a  special  method. 

1.3.1  The  Gauss-Newton  Method 

Vanilla  gradient  descent  is  the  simplest  and  most  intuitive  technique  to  find  minima  in  a  function, 
but  it  suffers  from  various  convergence  problems.  This  situation  can  be  improved  upon  by  using 
curvature  as  well  as  gradient  information,  namely  second  derivatives.  One  way  to  do  this  is  to 
use  Newton’s  method  to  solve  the  equation  g(x)  =  J(x)Tf(x)  =  0. 

From  (1.3),  the  Newton  equation  becomes  G^pj.  =  -gfc,  or: 

(JfcJfc  +  Qfc)Pfc*“-Jfcffc  (1-4) 

where  a  quantity  subscribed  by  k  denotes  that  quantity  evaluated  at  xfe,  the  current  estimate  of 
the  solution.  Let  pN  denote  the  solution  of  (1.4),  i.e.,  the  Newton  direction. 

If  HffcH  tends  to  zero  as  x*,  approaches  the  solution,  matrix  Qk  also  tends  to  zero,  so  equa¬ 
tion  (1.4)  can  be  approximated  by  the  following  equation,  which  involves  only  the  first  derivatives 
of  f: 

3l3kPk  =  -Jlfk  (1.5) 

The  solution  of  (1.5)  is  a  solution  of  the  linear  least-squares  problem 

minimize  i  ||  Jfcp  +  f*,  ||^  (1.6) 

and  is  unique  if  3k  has  full  column  rank.  The  vector  pGN  that  solves  (1.6)  is  called  Gauss- 
Newton  direction,  and  the  method  in  which  this  vector  is  used  as  a  search  direction  is  known 
as  the  Gauss-Newton  method.  Early  implementations  of  the  Gauss-Newton  method  typically 
formed  the  explicit  matrix  JkJk  and  computed  pGN  by  solving  equation  (1.5).  The  disadvantage 
of  this  approach  is  that  the  condition  number  of  JkJk  is  the  square  of  that  of  Jk.  Consequently, 
unnecessary  error  may  occur  in  determining  the  search  direction. 

Ill-conditioning  is  a  common  feature  of  nonlinear  least-squares  problems  derived  from  pa¬ 
rameter  estimation  problems,  because  the  underlying  mathematical  model  is  often  ill-defined. 
Unnecessary  worsening  of  the  conditioning  can  be  avoided  by  solving  the  linear  least-squares 
problem  (1.6)  using  the  complete  orthogonal  factorization  or  the  singular-value  decomposition. 

In  the  rank-deficient  case,  any  implementation  that  uses  a  minimum-norm  solution  to  (1.6) 
must  include  a  strategy  for  estimating  the  rank  of  J*,. 
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1.3.2  The  Levenberg-Marquardt  Method 

It  can  be  seen  that  simple  gradient  descent  and  Gauss-Newton  iteration  are  complementary  in 
the  advantages  they  provide.  Levenberg  [19]  proposed  an  algorithm  based  on  this  observation, 
whose  update  rule  is  a  blend  of  the  above  mentioned  algorithms. 

The  Levenberg  search  direction  is  defined  as  the  solution  of  the  equations 

(J%Jk  +  \kI)Pk  =  -jZfk  (1.7) 

where  \k  is  a  non-negative  scalar. 

If  the  error  reduces  following  an  update,  it  implies  that  the  quadratic  assumption  on  F(x)  is 
working  and  the  scalar  \k  con  be  reduced  (usually  by  a  factor  of  10)  to  reduce  the  influence  of 
gradient  descent.  On  the  other  hand,  if  the  error  goes  up,  it  is  necessary  to  follow  the  gradient 
more;  thus,  \k  is  increased  by  the  same  factor. 

The  above  algorithm  has  the  disadvantage  that  if  the  value  of  \k  is  large,  the  calculated 
Hessian  matrix  is  not  used  at  all.  We  can  derive  some  advantage  out  of  the  second  derivative  even 
in  such  cases  by  scaling  each  component  of  the  gradient  according  to  the  curvature.  This  should 
result  in  larger  movement  along  the  directions  where  the  gradient  is  smaller.  This  crucial  insight 
was  provide  by  Marquardt  [21].  He  replaced  the  identity  matrix  in  (1.7)  with  the  diagonal  of  the 
Hessian  resulting  in  the  Levenberg-Marquardt  update  rule 

(JfcJfc  +  Afc  diag(Gfc))pfc  =  -Jfc  ffc  (1.8) 

Since  the  Hessian  is  proportional  to  the  curvature  of  F(x),  (1.8)  implies  a  large  step  in  the  direction 
with  low  curvature  and  a  small  step  in  the  direction  with  high  curvature. 

It  is  to  be  noted  that  while  the  LM  method  is  in  no  way  optimal  but  is  just  a  heuristic,  it 
works  extremely  well  in  practice.  The  only  flaw  is  its  need  for  matrix  inversion  as  part  of  the 
update.  Even  though  the  inverse  is  usually  implemented  using  clever  pseudo-inverse  methods 
such  as  singular  value  decomposition,  the  cost  of  the  update  becomes  prohibitive  after  the  model 
size  increases  to  a  few  thousand  parameters.  For  moderately  sized  models  (of  a  few  hundred 
parameters)  however,  this  method  is  much  faster  than  say,  vanilla  gradient  descent. 

Historically,  the  LM  algorithm  was  presented  by  Marquardt  as  given  above,  where  the  pa¬ 
rameter,  A,  was  manipulated  directly  to  find  the  minimum.  Subsequently,  a  trust-region  approach 
to  the  algorithm  has  gained  ground.  The  idea  of  the  model  trust-region  approach  is  to  accept 
the  minimum  of  the  quadratic  model  only  as  long  as  the  quadratic  model  adequately  reflects  the 
behavior  of  F.  Usually,  the  decision  as  to  whether  the  model  is  acceptable  is  based  on  the  norm 
of  the  computed  search  direction. 

A  unit  step  is  always  taken  along  pfc  in  (1.7),  i.e. ,  xfc+i  is  given  by  xfc  +  pfc,  where  pk  is  the 
solution  of  the  constrained  subproblem 

minimize  \  || Jfep  +  U  |||  (1.9) 

peK™  z 

subject  to  ||p||2  <  A  (1-10) 

for  some  A  >  0.  It  can  be  shown  that  the  solution  of  the  equations  (1.7)  solves  the  subproblem  (1.9) 
if  either  A  =  0  and  ||p||2  <  A,  or  A  >  0  and  ||p||2  =  A. 

Thus,  if  A  is  large  enough,  the  solution  of  (1.9)  is  simply  the  Newton  direction  (i.e.,  the 
solution  of  (1.7)  with  A  =  o).  Otherwise,  the  restriction  on  the  norm  will  apply,  and  ||p||2  =  A.  The 
search  direction  is  typically  found  by  solving  (1.9)  for  trial  values  of  A  and  evaluating  F  at  the 
resulting  trial  points.  A  vector  p  such  that  F(xfe+p)  is  sufficiently  less  than  Fk  must  exists  for  small 
enough  A  since  the  second-order  term  of  the  model  function  may  be  made  small  compared  to  the 
first-order  term.  As  A  ->  0,  ||p||2  0  and  p  becomes  parallel  to  the  steepest-descent  direction. 

1.4  Membrane  Inverse  Analysis 

A  membrane  is  a  thin  structural  element,  i.e.,  an  element  whose  size  in  one  direction  (the  thickness, 
z  or  3  in  the  following)  is  extremely  small  compared  with  the  other  two  ( x  and  y  or  1  and  2  in 
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Figure  1.1:  Membrane  model. 


the  following),  see  Figure  1.1.  The  transverse  gradients  of  strain  and  stress  components  along  the 
reference  surface  of  the  membrane  are  negligible.  As  a  consequence,  a  membrane  cannot  withstand 
transverse  loads  unless  some  in-plane  prestress  is  present. 

1.4.1  Kinematics 

Assuming,  from  the  point  of  view  of  kinematics,  that  in-plane  strains  are  uniform  throughout 
the  thickness  of  the  membrane,  strain  measurements  on  just  one  side  of  the  membrane  would 
be  sufficient  to  characterize  the  in-plane  strain  field  (on  the  contrary,  in  the  case  of  a  shell,  the 
average  of  the  measurements  on  both  surfaces  would  be  needed  to  eliminate  the  effect  of  bending) . 

The  membrane  strains,  i.e. ,  the  in-plane  components  of  the  Green-Lagrange  strain  tensor, 
defined  as: 


l/*  +  2  ( 

+  V%  +  w/x) 

(1.11a) 

'/»  +  \  ( 

u/y+v/y+w/y) 

(1.11b) 

^/y  ^  /  x 

“I”  ^/x^/y  V/xV/y  “h  W / x'm /yi 

(1.11c) 

can  be  collected  in  a  vector  e: 


or 


(1.12) 


£ij  ~ 


1 

2 


+  Uj/i  +  u  £ 


(1.13) 


where  (*)/(*)  indicates  the  derivative  of  (A)  with  respect  to  (4),  vector  u  =  {u;v;w}  =  {mi;«2;«3} 
collects  the  displacement  components  in  the  plane  of  the  membrane,  u  =  ui  and  v  =  «2  and  the 
one  along  the  transverse  direction,  w  =  u3. 


1.4.2  Cost  Function 

The  inverse  kinematics  problem  can  be  formulated  by  defining  an  appropriate  cost  function  of 
the  error  e  between  the  measured  and  the  configuration-dependent  strains. 

Consider  a  set  of  strain  measurements  e£m),  and  7^,  e.g.  DIC,  namely: 

(1.14) 
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that  correspond  to  the  configuration-dependent  strains  defined  earlier. 
The  error  e  is 


e  =  e(Grad(u))  -  (1.15) 

with  Grad(*)  =  {(*)/*;  (*)/„}. 

The  following  cost  function  is  considered: 

<t>(u)  =  f>e(u)  +  fc<t>u  (1.16) 

where: 

•  $e(u)  is  a  quadratic  function  of  the  error  e, 

<t>e(Grad(u))  —  -  /  eTDe  d A  (1.17) 

2  J  A 

with  D  an  arbitrary  positive  definite  weighting  matrix;  e.g.,  but  not  necessarily,  the  plane 
stress  constitutive  properties  matrix,  which  for  isotropic  materials  is 

E 


0 

0 

(1  -  v)/2 


(1.18) 


•  <t>u  is  a  regularization  contribution  in  the  derivatives  of  w, 

't’ufu)  =  —  f  (Grad(uj)  —  Grad(roref))TT(Grad(uj)  —  Grad(roref))  dT  (1.19) 

2  J  A 

with  wle{  a  reference  transverse  displacement,  defined  in  Section  1.4.3,  and  the  weighting 
matrix 


T  = 


(1.20) 


defined  in  analogy  with  the  strain  energy  contribution  associated  with  pretension:  Tx  >  0, 
Ty  >  0,  and  A /TxTy  >  \Txy  \  >  0  such  that  T  >  0  (positive  definite); 


•  k  is  a  parameter  that  restores  dimensional  consistency  and  weighs  the  regularization  con¬ 
tribution. 


1.4.3  Regularization 

The  regularization  contribution  is  defined  in  such  a  manner  that  it  naturally  vanishes  at  con¬ 
vergence,  by  properly  crafting  the  reference  displacement,  wre{.  Such  correction  is  needed  to  add 
a  positive  definite  quadratic  contribution  to  the  cost  function,  and  thus  make  it  convex  on  the 
entire  domain.  In  fact,  the  minimization  of  <J>e(u)  with  respect  to  the  actual  displacement  field 
u  requires  its  partial  derivatives  with  respect  to  each  of  the  components  u,  v,  and  w  to  vanish. 
Clearly,  as  a  consequence  of  the  strain  definitions  of  Eq.  1.11,  <E>e(Grad(u))  is  not  a  convex  function 
of  w  when  Grad(u>)  =  0,  i.e.,  when  the  membrane  is  parallel  to  the  reference  plane,  and  the  problem 
is  ill-posed  in  the  vicinity  of  such  condition. 

This  approach  resembles  the  so-called  damped  least  squares,  also  known  as  the  Levenberg- 
Marquardt  algorithm  [19,  21],  summarized  in  Section  1.3.2. 

The  inverse  formulation  does  not  need  elastic  or  inertial  material  properties.  The  reference 
transverse  displacement  wref  ^  0  is  needed  to  deflect  the  solution  towards  a  specific  direction, 
since  the  same  membrane  strain  pattern  is  obtained  with  ±w.  As  suggested  in  Alioli  et  al.  [4],  one 
should  choose  a  tentative  initial  value  for  turef:  a  convenient  choice  can  be  the  one  corresponding  to 
a  uniform  pressure  difference  applied  on  the  membrane,  or  in  any  case  a  prescribed  displacement 
that  qualitatively  resembles  the  expected  solution.  Subsequently,  the  reference  solution  is  updated 
by  interpolating  between  the  current  value  and  the  solution  at  the  current  step  i,  i.e.,  «/')  + 
Auidl,  such  that  at  convergence,  when  Au>W  =  0,  then  (Grad(w)  -  Grad(u>ref))  =  0.  Thus,  = 

(1  —  «)^ref  +  a(w^  +  Auid)),  where  0  <  a  <  1  is  a  relaxation  parameter  (a  =  1  implies  no  relaxation). 
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Figure  1.2:  Deformation  error  of  a  square  membrane  subjected  to  uniform  pressure. 


1.4.4  Pressure  Field  Reconstruction 

Finally,  the  estimated  displacement  field  u  =  =  {«i;  112;  113}  can  be  used  to  estimate  the 

distributed  force  field  p  acting  on  the  membrane  [6]: 

f  Se  :  <j  dA  =  f  p  ■  Sv  d A  (1.21) 

J  A  J  A 

where  <r  =  <r(e)  is  the  known  stress  tensor,  expressed  as  a  function  of  the  reconstructed  strain 
tensor  e  =  e(u),  and  Sv  is  an  appropriate  vector  test  function ,  which  is  required  to  vanish  on  the 
boundary  8A  of  A.  Thus,  the  reconstructed  distributed  force  field  p  is  expected  not  to  be  accurate 
along  the  boundary  of  the  membrane  domain,  due  to  the  choice  of  the  boundary  conditions. 

1.4.5  Verification 

In  this  section,  a  simple  direct  ( “forward” )  membrane  problem  is  considered  in  order  to  demon¬ 
strate  how  the  inverse  formulation  lead  to  the  reconstruction  of  the  displacement  solution  consis¬ 
tent  with  the  measured  strain  data.  The  exact  strains  are  computed  from  the  direct  problem  and 
are  used  to  represent  the  measured  strain  data. 

Consider  a  direct  (forward)  problem  of  a  rectangular  membrane  bent  under  the  action  of  the 
distributed  transverse  loading  (hydrostatic  pressure).  The  proposed  inverse  finite  element  analysis 
formulation  is  verified  in  Alioli  et  al.  [4,  5]  performing  a  direct  analysis,  obtained  with  the  same 
finite  element  model  (an  edge-clamped  square  membrane  without  prestrain  and  subjected  to 
uniform  pressure  on  one  side),  but  with  a  finer  mesh,  whose  results  are  re-sampled  to  provide  the 
membrane  strain  values  at  the  points  required  by  the  inverse  analysis.  The  re-sampled  transverse 
displacements  are  used  to  verify  the  quality  of  the  inverse  analysis.  The  numerical  solution, 
obtained  using  a  fine  mesh  (20x20  pairs  of  elements),  has  been  re-sampled  using  a  much  coarser 
mesh  (10x10).  A  third,  intermediate  mesh  (12x12)  has  been  used  for  the  IFEM  procedure. 

Figure  1.2  shows  the  transverse  displacement  error,  computed  as  weTT  =  |w(fem)  —  w|/max  |w(fem)|, 
where  w(£em 1  is  the  transverse  displacement  obtained  by  the  reference  finite  element  analysis, 
whereas  w  is  the  transverse  displacement  reconstructed  via  IFEM. 

The  estimated  displacement  field  is  subsequently  used  to  estimate  the  pressure  acting  on  the 
membrane:  Table  1.1  compares  the  resultant  (normal)  force  error  obtained  by  using  the  FEM 
pressure  distribution  (which  is  calculated  from  the  displacements  of  the  reference  FEM  analysis) 
and  the  one  obtained  by  using  the  pressure  distribution  computed  from  the  displacements  esti¬ 
mated  via  IFEM  analysis,  with  respect  to  the  resultant  (normal)  force  obtained  considering  the 
ideal  uniform  pressure  distribution,  i.e.  pressure  times  the  membrane  surface  area.  To  compute 
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the  resultant  forces,  the  near-boundary  region  was  omitted  from  calculation,  for  the  motivation 
just  mentioned  above. 


Reference  FEM  solution  1.64  % 
Inverse  FEM  analysis  5.25  % 


Table  1.1:  Resultant  normal  force  error  (not  considering  the  near-boundary  region),  comparison  of  FEM 
and  IFEM  analysis  with  respect  to  the  nominal  value. 

The  results,  in  terms  of  transverse  displacement,  are  satisfactory.  Furthermore,  as  shown 
in  Alioli  et  al.  [6],  the  IFEM  and  the  FEM  pressure  distributions  are  quite  similar  on  the  entire 
domain,  although  both  of  them  are  not  significantly  flat.  This  is  expected,  due  to  the  choice  of 
the  boundary  conditions.  Nevertheless,  the  average  values  of  the  pressure  distributions  calculated 
using  the  reconstructed  displacements  and  those  from  the  FEM  analysis  do  seem  to  be  accurate, 
as  seen  from  the  resultant  force  errors  (Tab.  1.1). 

A  typical  inverse  solution  like  the  ones  presented  above,  e.g.,  with  a  12  x  12  mesh  and  a  residual 
norm  tolerance  of  10-5,  requires  5  to  8  iterations  for  each  load  step.  Each  iteration  requires  about 
27.5  ms  on  an  off-the-shelf  PC  (in  the  present  case,  an  Intel  Core  i7-2620M  with  CPU  at  2.70  GHz). 
The  property  of  computational  efficiency  is  of  utmost  importance  since  the  long  term  objective 
is  the  real-time  implementation  of  the  procedure. 

1.4.6  Problem  Well-Posedness 

As  we  shall  see,  the  membrane  inverse  problem  may  be  ill-posed,  and  can  have  multiple  solutions. 

The  geometry  of  a  two-dimensional  surface  is  characterized  by  two  symmetric  tensors,  called 
the  first  and  second  fundamental  forms  of  a  surface.  See  Do  Carmo  [13]  for  details. 

Let  x(£“),  with  a  =  {1,  2},  be  the  parametric  equations  of  the  surface.  Let  also  ga  =  x:a  =  ^  and 
g“  be  the  covariant  and  contravariant  surface  base  vectors,  respectively,  defined  in  such  a  way 
that  ga  •  g'3  =  <5^,  the  Kronecker  delta,  with  a,  /3  =  {1,  2}. 

The  first  fundamental  form  a  of  the  surface  is  nothing  but  its  metric  tensor,  defined  in  such  a 
way  that  x,a  =a  ga.  It  can  be  computed  as  a  =  g“  ig>  ga  •  g^  <s>  g'3  ■  The  second  fundamental  form  b 
of  the  surface  describes  the  rate  of  change  of  the  surface  normal  n,  viz.:  n,a  =  — b  •  gQ. 

The  fundamental  theorem  of  the  theory  of  surfaces  states  that  the  first  and  second  fundamental 
forms  of  a  surface  determine  its  shape  up  to  its  position  in  space  (for  a  self-contained,  and 
essentially  elementary  proof,  see  Ciarlet  and  Larsonneur  [12]).  Thus,  this  theorem  guarantees  that, 
by  knowing  both  the  first  and  the  second  fundamental  form,  the  inverse  problem  of  reconstructing 
the  surface  is  well-posed. 

The  ill-posedness  of  the  inverse  problem  comes  from  the  fact  that,  by  measuring  only  the 
membranal  strain  tensor  e  one  actually  accounts  only  for  the  first  fundamental  form,  and  not 
for  the  second.  As  a  matter  of  fact,  the  membrane  strain  tensor  e,  Eq.  (1.12),  is  nothing  but  the 
difference  between  the  metric  tensor  a'  in  the  deformed  configuration  x'  =  x  +  u,  computed  with 
respect  to  the  reference  configuration,  and  the  metric  tensor  in  the  reference  configuration  a,  viz.: 
s  —  l(a'  —  a).  The  difference  between  the  second  fundamental  form  in  deformed  and  reference 
configuration,  i.e.,  k  =  b'  —  b,  is,  instead,  a  suitable  measure  of  the  flexure  strain  of  a  thin  shell. 
This  means  that  the  inverse  problem  would  be  well-posed  for  a  shell  model,  provided  both  the 
membrane  and  flexure  strain  are  measured. 

As  an  example  of  an  inverse  problem  that  is  not  solvable,  consider  a  rectangular  membrane 
made  of  isotropic  material,  with  Poisson  coefficient  v  =  0,  and  subjected  to  a  cylindrical  bending 
with  a  constant  deformation  eh  =  const  and  all  the  other  deformation  components  equal  to  zero, 
i.e.  £22  =  ei2  =  0.  It  is  trivial  to  verify  that  this  problem  has  not  a  unique  solution.  For  example, 
both  the  deformed  configurations  of  Figure  1.3  are  among  the  possible  solutions  of  this  problem. 

At  this  point  one  could  ask  himself  whether  the  inverse  procedure  described  in  this  report 
leads  to  meaningful  results,  if  any.  And,  indeed,  the  procedure  fails,  as  it  should,  when  applied 
to  the  ill-posed  cylindrical  bending  problem  described  above.  Nonetheless,  it  appears  to  work 
reasonably  well  for  the  test  cases  that  were  defined  for  this  activity. 
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Figure  1.3:  Two  possible  solutions  for  a  cylindrical  bending  inverse  problem. 


An  insight  into  this  apparent  paradox  comes  from  the  fact  that  the  covariant  base  vectors  gQ 
can  be  used  to  compute  the  surface  normal  as  n  =  Then,  since  g„  ,g7  =  g arfp,  there  must 

be  some  compatibility  equations  between  the  first  and  the  second  fundamental  forms.  These  equa¬ 
tions,  known  as  Codazzi-Gauss  equations,  can  be  written  in  many  ways.  One  of  them,  reported 
below,  is 


611.2  —  612,1  =  0  (1.22a) 

621.2  —  622,1  =  0  (1.22b) 

bub22-b212  =  Ka  (1.22c) 

where  K  =  det(b)  is  the  Gaussian  curvature  of  the  surface  and  a  =  det(a).  Thus,  what  makes  some 
problems  solvable  is  the  fact  that  the  deformed  configuration  has  a  non-null  Gaussian  curvature, 
K  =  det(b)  ^  0,  so  that  a  link  can  be  implicitly  established  between  the  first  and  the  second 
fundamental  form  of  the  surface. 

From  a  physical  point  of  view,  what  is  really  important  is  that  the  deformed  configuration 
should  not  allow  an  additional  bending  deformation  of  the  surface  that  does  not  involve  an 
additional  membrane  deformation  as  well.  In  other  words,  the  problem  should  be  membrane- 
dominated,  in  the  sense  defined  by  Chapelle  and  Bathe  [11]. 

1.5  Direct  Analysis 

The  direct  analysis  is  performed  using  a  tightly  coupled  fluid-structure  co-simulation  in  which  the 
structural  problem  is  solved  using  the  free  general-purpose  multibody  dynamics  solver  MBDyn2 
[23]  and  the  fluid  problem  is  solved  using  a  dedicated  solver  based  on  FEniCS  [2,  3],  where  systems 
of  Partial  Differential  Equations  (PDE)  and  corresponding  discretization  and  iteration  strategies 
can  be  defined  in  terms  of  a  few  high-level  Python  statements  which  inherit  the  mathematical 
structure  of  the  problem,  and  from  which  low  level  code  is  automatically  generated.  The  fluid 
dynamics  code  is  based  on  a  stabilized  finite  element  approximation  of  the  unsteady  Navier-Stokes 
equations  (often  referred  to  in  the  literature  as  G2  method  [16]). 

The  multibody  solver  is  coupled  with  the  external  fluid  dynamic  code  by  means  of  a  general- 
purpose,  meshless  boundary  interfacing  approach  based  on  Moving  Least  Squares  with  Radial 
Basis  Function,  as  presented  in  Quaranta  et  al.  [24] .  This  technique  allows  to  compute  a  sufficiently 
regular  and  accurate  approximation  of  the  field  of  the  structural  displacements  and  velocities  at 
the  aerodynamic  interface  nodes,  based  on  a  set  of  structural  nodes  that  is  in  general  irregularly 
distributed  in  the  neighborhood  of  the  interface. 

The  membrane  element,  implemented  in  MBDyn  as  shown  in  Masarati  et  al.  [22],  is  formulated 
as  a  four-node  isoparametric  element  based  on  second  Piola-Kirchhoff  type  membranal  resultants. 
The  classical  Enhanced  Assumed  Strains  (EAS)  method  [27]  is  exploited  to  improve  the  response 

2http: //www.mbdyn. org/. 
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of  the  element:  seven  additional  variables  for  each  membrane  element  are  added  to  the  strain 
vector  (see  for  example  Andelfmger  and  Ramm  [8]  for  details). 

The  stress  tensor,  reorganized  in  form  of  a  vector,  can  be  expressed  as  a  function  of  the  strain 
tensor,  reorganized  in  the  same  manner,  using  the  constitutive  law  of  the  membrane  element,  e.g., 
Eq.  (1.18): 


(1.23) 


In  case  of  homogeneous  constitutive  properties,  the  forces  per  unit  span  are  readily  obtained  by 
multiplying  the  stresses  by  the  thickness  h  of  the  membrane  (otherwise,  thickness-wise  integra¬ 
tion  is  required).  Generically  anisotropic  constitutive  properties  can  be  defined,  with  matrix  D 
symmetric,  positive  definite  but  otherwise  arbitrarily  set  by  the  user. 
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This  section  presents  the  validation  of  the  direct  and  inverse  analysis  by  comparing  the  direct 
prediction  and  the  reconstructed  deformed  shape  with  experimental  data  for  prestressed  rectan¬ 
gular  membrane  wings  subjected  to  hydrostatic  pressure  loads  (Section  2.2)  and  in  steady  level 
flight  (Section  2.3).  Since  the  methods  described  above  are  being  evaluated  for  their  capability  of 
estimating  an  actual  pressure  distribution,  it  is  first  necessary  to  know  what  the  applied  pressure 
is,  in  order  to  have  a  basis  for  comparison.  For  this  reason,  two  experiments  were  conducted  in 
this  work  to  provide  different  loading  scenarios  for  the  estimation  routine. 

The  first  scenario  was  a  hydrostatic  pressure  case,  where  a  pre-tensioned  membrane  was  sub¬ 
jected  to  a  constant  known  pressure.  DIC  measurements  of  strain  and  deformation  were  taken  of 
the  deformed  membrane. 

In  the  second  scenario,  a  membrane  wing  was  placed  in  a  low  speed  wind  tunnel,  and  wind 
speed  and  angle  of  attack  were  varied.  Aerodynamic  loads  generated  by  the  wing  were  measured, 
and  DIC  measurements  of  the  membrane  deformation  were  taken. 

Correlation  is  sought  with  respect  to  experimental  results  obtained  in  test  campaigns  per¬ 
formed  at  Oregon  State  University,  where  elastic  deformations  and  strains  were  measured  using 
DIC  [10]. 

Figure  2.1  summarizes  the  verification  procedure:  the  strain  measurements  are  first  re-sampled 
onto  the  numerical  grid  that  is  subsequently  used  for  IFEM  analysis  by  means  of  the  previously 
discussed  MLS  procedure  using  radial  basis  functions  (initially  developed  for  field  interpolation 
at  the  interface  between  fluid  and  structure,  see  Quaranta  et  al.  [24]  for  further  details).  The 
re-sampled  measurements  are  used  as  inputs  for  the  IFEM  analysis.  The  re-sampled  (transverse) 
displacements  are  used  to  evaluate  the  quality  of  the  reconstructed  displacements  via  IFEM  and 
those  predicted  using  the  tightly  coupled  fluid-structure  co-simulation. 


Figure  2.1:  Scheme  of  the  verification  procedure. 

Furthermore,  total  force  measurements  are  used  to  evaluate  the  total  force  reconstructed  by 
direct  analysis  and  IFEM,  and  pressure  distributions  determined  by  direct  analysis  are  used  to 
evaluate  the  pressure  distribution  reconstructed  by  IFEM. 
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2.1  Experimental  Data  Re-Sampling 

Measurements  provided  by  DIC  [10]  include:  (i)  the  reference  location  in  space  of  an  arbitrary  set 
of  points  on  the  surface,  chosen  by  the  DIC  algorithm  when  the  measurement  system  is  activated, 
(ii)  the  displacements  of  the  corresponding  points  in  the  current  sample,  (iii)  an  estimate  of  the 
in-plane  strains. 

Data  preparation,  for  both  the  measured  strains  used  as  inputs  and  the  measured  displace¬ 
ments  used  for  correlation,  requires  re-sampling  of  unstructured  measured  fields  onto  the  grid  that 
is  subsequently  used  for  IFEM  analysis.  This  is  done  to  reduce  the  size  of  the  IFEM  problem,  to 
avoid  distortions  in  the  mesh  of  the  IFEM  model,  and  to  obtain  an  initial  spatial  filtering  of  the 
measurements . 

For  this  purpose,  a  meshless  mapping  procedure  originally  developed  for  fluid-structure  cou¬ 
pling  is  used.  The  mapping  [24]  produces  a  linear  interpolation  operator,  W,  from  the  measurement 
domain,  (-)m,  to  the  virtual  sensing  domain,  (•)„,  namely  x„  =  H-Xm-  Operator  H  is  computed 
based  on  the  initial  positions  of  both  domains;  from  that  point  on,  it  is  used  to  map  an  arbitrary 
configuration  of  the  measure  domain  onto  the  virtual  sensing  domain. 

The  participation  of  each  component  of  a  measure  point’s  position  to  the  mapping  of  the 
corresponding  component  of  a  virtual  point  is  the  same,  i.e. ,  the  mapping  is  isotropic.  As  a 
consequence,  any  scalar  field,  as  well  as  each  component  of  any  vector  field,  can  be  mapped 
separately  using  a  subset  of  matrix  H,  obtained  for  example  by  extracting  every  one  out  of  three 
columns  and  rows  of  matrix  H, 


H  =  W(l:3:end,  l:3:end).  (2.1) 

The  component-by-component  re-sampled  strain  measurements,  eiv  =  Heim,  i  =  11,22,12,  are 
used  as  inputs  for  the  IFEM  procedure,  whereas  the  re-sampled  displacements,  u„  =  4tum,  are 
used  to  evaluate  the  quality  of  the  reconstructed  displacements. 

When  a  sufficiently  large  number  of  measurement  points  is  required  to  interpolate  the  position 
of  a  virtual  sensing  point,  as  occurs  in  the  present  case,  the  procedure  also  produces  a  smoothing 
of  the  input  data,  acting  as  a  spatial  filter. 

Thanks  to  the  compact  support  used  for  the  interpolation  [24],  the  mapping  matrix  H  is 
usually  quite  sparse:  Fig.  2.2  shows  the  shape  and  fill-in  (of  the  order  of  0.05  %  of  non-zeroes)  for 
the  4x8  (45  nodes),  8  x  16  (153  nodes)  and  16  x  32  (561  nodes)  membrane  meshes  mapped  from 
6416  DIC  points.  Matrix  H  can  be  (and  it  is,  indeed)  stored  and  handled  exploiting  such  sparsity, 
thus  drastically  reducing  the  computational  cost  associated  with  field  mapping  (see  Alioli  et  al. 
[4]  for  details). 

2.2  Hydrostatic  Pressure  Test 

2.2.1  Test  Article 

The  experiments  in  Carpenter  and  Albertani  [10]  refer  to  a  rectangular  edge-clamped  membrane 
wing,  whose  dimensions  are  140  x  75  x  0.14  mm.  The  wing  was  constructed  of  a  pretensioned  rubber 
latex  membrane.  The  material  properties  are  reported  in  Table  2.1.  The  membrane,  prestrained 
by  a  9%  isotropic  membrane  strain  ( eXQ  =  eyo  =  0.09,  ^xyo  =  o),  was  subjected  to  hydrostatic 
pressure  difference  between  the  lower  and  the  upper  surface  ranging  from  100  Pa  to  500  Pa  in 
steps  of  100  Pa. 


Tensile  modulus  E 

Poisson’s  modulus  v 

Density  p 

latex  rubber  1.8  MPa 

0.4 

1350.  kg/m3 

Table  2.1:  Membrane  material  properties. 
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Figure  2.2:  Shape  and  fill-in  of  mapping  operator  H. 
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■  MBDyn 
•  mapped  DIC 


(a)  300  Pa. 


■  MBDyn 
•  mapped  DIC 


(b)  500  Pa. 


Figure  2.3:  Numerical/experimental  correlation  of  hydrostatic  pressure  problem. 


2.2.2  Direct  Survey 

Figures  2.3(a)  and  2.3(b)  compare  the  numerical  results,  for  the  problem  with  300  Pa  and  500  Pa 
of  pressure  difference,  with  the  experimental  ones  re-sampled  on  the  numerical  mesh  using  the 
previously  discussed  moving  least  squares  (MLS)  procedure.  The  same  domain  mapping  algo¬ 
rithm  is  used  to  exchange  motion  and  loads  at  the  nodes  between  MBDyn  and  the  fluid  solver 
implemented  in  FEniCS  during  the  coupled  fluid-structure  solution  when  the  interface  nodes  of 
the  structure  and  fluid  domains  do  not  match. 

The  structural  grid,  implemented  within  the  multibody  simulation  environment  provided  by 
MBDyn,  consists  of  8x16  four-node  membrane  elements,  involving  153  structural  nodes  (and 
153  rigid  body  elements  when  a  dynamic  model  needs  to  be  used),  as  shown  in  Figure  2.4(a). 
Although  not  involved  in  the  presented  test  cases,  the  mass  lumped  in  each  node  is  computed 
from  the  latex  rubber  sheet  portion  associated  with  the  node,  which  is  uniformly  distributed,  see 
Table  2.1. 

A  comparison  of  the  direct  solution  for  smaller  and  larger  values  of  the  prestrain  is  performed 
to  study  how  much  the  problem  is  dependent  on  the  value  of  the  prestrain  [4,  6].  In  the  experiments 
described  in  Carpenter  and  Albertani  [10],  the  membrane  prestrain  was  introduced  as  accurately 
as  possible  at  the  nominal  level,  but  could  not  be  checked  afterwards.  The  facts  that  the  numerical 
solutions  with  nominal  prestrain  present  a  very  good  correlation  with  the  experiments,  and  that 
the  solution  is  very  sensitive  to  the  amount  of  prestrain,  indicate  that  the  actual  prestrain  in  the 
experiments  was  in  accordance  with  the  expected  value.  The  inverse  solution,  instead,  is  insensitive 
to  the  amount  of  prestrain  (matrix  T  of  Eq.  (1.20))  that  is  used  to  stabilize  the  solution  process. 


(a)  Membrane  model.  (b)  Composite  model. 

Figure  2.4:  Multibody  membrane  models. 
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is  average  model 

■  average  model 

•  fiber  model 

•  fiber  model 

(a)  100  Pa.  (b)  200  Pa. 

Figure  2.5:  Comparison  of  the  “fiber”  and  “average”  model,  hydrostatic  pressure  problem. 


Non-Isotropic  Membrane 

A  hydrostatic  membrane  pressure  test  was  also  conducted  to  characterize  the  behavior  of  a 
non-isotropic  membrane  under  constant  uniform  pressure  distribution.  A  study  case  of  non- 
homogeneous  constitutive  properties  of  the  model  was  investigated  using  the  100%  fiber  modulus 
from  experimental  data  [32].  The  membrane  dimensions,  loading,  and  boundary  conditions  are 
the  same  as  in  the  direct  analysis  case  discussed  above.  No  prestrain  was  introduced  herein. 

A  composite  membrane  was  considered:  seven  fibers,  which  are  modeled  as  “rod”  elements, 
are  oriented  parallel  to  the  spanwise  direction  (i.e. ,  the  y-axis),  as  shown  in  Figure  2.4(b).  The 
composite  is  made  of  spandex  for  the  fibers  and  silicone  for  the  matrix.  Table  2.2  presents  the 
properties  of  each  material.  The  elastic  modulus  at  100%  strain  was  matched  with  experimental 
results  of  tensile  tests  [32].  The  fiber  diameter  is  equal  to  the  thickness  of  the  membrane. 

A  corresponding  model  of  the  specimen  with  averaged  membrane/fiber  properties  was  also 
investigated,  by  defining  orthotropic  constitutive  properties  of  the  membrane  finite  element. 


Tensile  modulus  E 

Poisson’s  modulus  v 

silicone  matrix 

0.379  MPa 

0.4 

spandex  fiber 

2.2  MPa 

0.4 

Table  2.2:  Matrix  and  fiber  material  properties. 

Figures  2.5(a)  and  2.5(b)  respectively  compare,  for  the  problem  with  100  Pa  and  200  Pa  of 
pressure  difference,  the  deformed  shapes  obtained  using  a  model  of  the  specimen  that  explic¬ 
itly  models  the  fibers  (the  “fiber”  model)  with  those  resulting  from  a  corresponding  model  in 
which  averaged  membrane/fiber  properties  were  used  (the  “average”  model).  As  shown,  the  two 
models  give  consistent  numerical  results  for  hydrostatic  pressure  tests  under  controlled  boundary 
conditions. 

2.2.3  Inverse  Survey 

The  deformed  configuration  of  the  previously  investigated  rectangular  edge-clamped  membrane 
is  determined  herein  using  the  strain  measurements  derived  from  DIC,  re-sampled  onto  the  nu¬ 
merical  mesh  using  the  previously  discussed  MLS  procedure.  The  membrane  dimensions,  loading, 
boundary  conditions,  and  material  properties  are  the  same  as  in  the  direct  analysis  case  discussed 
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above,  see  Table  2.1.  A  triangular  mesh  consisting  of  8x16  elements  has  been  used  for  the  IFEM 
procedure  (Figure  2.6):  the  strain  measurements  from  DIC  are  re-mapped  (and  smoothed  out) 
onto  this  virtual  strain  measurement  grid,  and  used  as  inputs  for  the  IFEM  analysis. 


Figure  2.6:  Membrane  inverse  analysis  mesh. 

In  Figs.  2.7(a)  and  2.7(b)  the  deformation  shape  corresponding  to  the  inverse  FEM  analysis 
for  the  problems  with  300  Pa  and  500  Pa  of  pressure  difference  are  presented,  along  with  the 
experimental  data  re-sampled  on  the  numerical  mesh,  to  validate  the  deformations  predicted  by 
the  IFEM  analysis. 


DIC 

DIC 

•  mapped  DIC 

•  mapped  DIC 

♦  iFEM 

♦  iFEM 

(a)  300  Pa.  (b)  500  Pa. 


Figure  2.7:  IFEM/DIC  correlation  of  hydrostatic  pressure  problem. 

In  order  to  evaluate  the  quality  of  the  numerical  solution,  a  study  of  the  sensitivity  of  the 
deformation  to  the  refinement  of  the  mesh  used  for  the  IFEM  analysis  is  performed  in  Alioli  et  al. 
[4].  In  this  case,  the  transverse  displacement  error  is  computed  as  welI  =  -  ru)/max(u!lml), 

where  w is  the  measured  transverse  displacements,  whereas  w  is  the  transverse  displacement 
reconstructed  via  IFEM.  As  shown  in  Figure  2.8,  the  problem  with  500  Pa  of  pressure  difference 
appears  to  be  essentially  at  convergence  even  with  a  mesh  consisting  of  8x16  elements,  and  the 
experimental  displacements  and  those  reconstructed  using  IFEM  are  in  good  agreement  on  the 
entire  domain,  including  the  maximum  values. 

The  internal  stresses  for  one  particular  hydrostatic  pressure  test  at  500  Pa  are  shown  in  Fig¬ 
ure  2.9.  The  stress  components  <jxx  and  ayy  grow  from  the  initial  value,  about  0.238  MPa,  to  a 
maximum  value  of  the  order  of  0.289  MPa  for  axx  and  0.276  MPa  for  ayyi  respectively,  in  the  central 
portion  of  the  membrane,  where  the  extension  is  maximal.  The  cross  term  axy  remains  negligible. 
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Chapter  2.  Results 


2.2.  Hydrostatic  Pressure  Test 


Figure  2.8:  Transverse  displacement  error  wer r  of  a  rectangular  membrane  subjected  to  500  Pa  of  pressure 
difference  (left:  16x32;  center:  8x16;  right:  4x8). 
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Figure  2.9:  Membrane  stress  distribution  (MPa)  of  a  rectangular  membrane  subjected  to  500  Pa  of 
pressure  difference. 
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Chapter  2.  Results 


2.3.  Wind  Tunnel  Tests 


Such  a  change,  of  the  order  of  20%,  clearly  shows  that  the  approximation  of  constant  membrane 
stress,  which  is  at  the  roots  of  linearized  membrane  model,  is  not  applicable  to  problems  of  this 
type. 

The  pressure  distribution  for  one  particular  hydrostatic  pressure  test  at  300  Pa  is  shown  in 
Figure  2.10,  in  order  to  compare  a  known  input  data,  in  this  case  a  static  pressure,  with  the 
estimated  load.  In  fact,  with  a  known  input  to  the  system,  the  output  from  the  loads  estimation 
procedure  could  be  directly  evaluated  for  its  accuracy.  The  Figure  on  the  left  shows  the  ideal 
hydrostatic  pressure  applied  to  the  membrane.  The  center  Figure  shows  the  pressure  distribution 
calculated  from  full  field  DIC  measurements  remapped  onto  the  numerical  mesh  grid  as  required 
by  the  IFEM  procedure.  The  Figure  on  the  right  shows  the  estimated  pressure  distribution  from 
the  estimated  displacements  via  IFEM  analysis.  The  results  from  the  hydrostatic  pressure  test 


Figure  2.10:  Pressure  distribution  (MPa)  of  a  rectangular  membrane  subjected  to  300  Pa  of  pressure 
difference  (left:  ideal  pressure  distribution;  center:  re-mapped  DIC  pressure  distribution;  right:  IFEM 
pressure  distribution). 


show  favorable  results:  the  average  hydrostatic  pressure  estimates  are  reasonably  close  to  the 
actual  applied  hydrostatic  load,  and  the  error  between  the  resultant  (normal)  force  from  the 
estimated  pressure  distributions  and  the  one  from  the  ideal  hydrostatic  pressure  distribution  is 
relatively  small,  as  shown  in  Table  2.3  and  in  Alioli  et  al.  [7]. 


Re-mapped  DIC  solution  0.42  % 
Inverse  FEM  analysis  8.1  % 


Table  2.3:  Resultant  normal  force  error  (not  considering  the  near-boundary  region),  comparison  of  re¬ 
mapped  DIC  and  IFEM  analysis  with  respect  to  the  nominal  value. 


2.3  Wind  Tunnel  Tests 

2.3.1  Test  Article 

The  experiments  conducted  at  OSU  also  involved  wind  tunnel  tests  of  various  2:1  aspect  ra¬ 
tio,  rectangular,  perimeter  reinforced  membrane  wings.  The  membrane  dimensions  and  the  ma¬ 
terial  properties  are  same  as  in  Section  2.2,  see  Table  2.1.  The  wing  was  constructed  of  two 
shaped  steel  frames,  sandwiching  a  pretensioned  rubber  latex  membrane.  The  overall  geometry 
was  140  mm  x  75  mm,  with  a  frame  width  and  thickness  of  5  mm  and  1  mm  respectively. 

The  wind  tunnel,  depicted  in  Figure  2.11(a)  and  2.11(b),  that  was  used  to  conduct  all  ex¬ 
periments  is  a  low-speed  wind  tunnel  located  at  Oregon  State  University,  Corvallis  (OR).  The 
wind  tunnel  had  a  closed  loop,  closed  test  section,  capable  of  speeds  from  1  to  18  meters  per 
second  (m/s)  and  with  a  1.3  x  1.5  m  test  section. 


Real-Time  Pressure  Distribution  Estimation  on  Wings  Via  Displacements  and  Strains 


20 


DISTRIBUTION  A.  Approved  for  public  release:  distribution  unlimited. 


Chapter  2.  Results 


2.3.  Wind  Tunnel  Tests 


The  parameters  considered  in  these  tests  were  AoA,  initial  prestrain  and  flow  velocity.  Wind 
tunnel  tests  are  run  for  three  different  feasible  MAV  flight  speeds  (12,  15  and  18  m/s),  at  three 
pre-stall  angles  of  attack  (3,  6  and  9  deg),  with  three  different  initial  prestrain  values  (2,  3.5  and 
5%).  The  maximum  Reynolds  number  is  67000.  At  each  flight  condition,  the  aerodynamic  loads 
are  measured  with  a  six  component  sting  balance.  At  the  same  time,  the  undeformed  wing  shape 
and  the  strain  field  are  measured  using  DIC. 

To  evaluate  the  validity  of  the  purposed  approach,  experimental  wind  tunnel  loads  and  DIC 
displacements  are  compared  to  those  obtained  with  the  purposed  model  under  varying  conditions 
of  flow  velocity,  AoA  and  initial  prestrain. 


(a)  Overall.  (b)  Zoom. 


Figure  2.11:  Wind  Tunnel  Apparatus. 


2.3.2  Direct  Survey 

Figures  from  2.12  to  2.15  compare  the  numerical  results,  with  the  experimental  ones  re-sampled 
on  the  numerical  mesh  using  the  previously  discussed  moving  least  squares  (MLS)  procedure. 
The  same  domain  mapping  algorithm  is  used  to  exchange  motion  and  loads  at  the  nodes  between 
MBDyn  and  the  fluid  solver  implemented  in  FEniCS  during  the  coupled  fluid-structure  solution. 
The  structural  grid  consists  of  8x16  four-node  membrane  elements,  involving  153  structural  nodes, 
and  thus  153  rigid  body  elements,  as  shown  in  Figure  2.4(a):  the  mass  lumped  in  each  node  is 
computed  from  the  latex  rubber  sheet  portion  associated  with  the  node,  which  is  uniformly 
distributed,  see  Table  2.1. 
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■  MBDyn 

■  MBDyn 

•  mapped  DIC 

•  mapped  DIC 

(a)  AoA  =  3  deg. 

Figure  2.12:  DIC/MBDyn  correlation  of  wind  tunnel  tests 


(b)  AoA  =  6  deg. 


(V  =  12  m/s  and  PS 


2%). 


■  MBDyn 

■  MBDyn 

•  mapped  DIC 

•  mapped  DIC 

(a)  AoA  =  3  deg.  (b)  AoA  =  6  deg. 


Figure  2.13:  DIC/MBDyn  correlation  of  wind  tunnel  tests  (V  =  15  m/s  and  PS  =  2%). 
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2.3.  Wind  Tunnel  Tests 


■  MBDyn 

■  MBDyn 

•  mapped  DIC 

•  mapped  DIC 

(a)  AoA  =  6  deg. 


(b)  AoA  =  9  deg. 


Figure  2.14:  DIC/MBDyn  correlation  of  wind  tunnel  tests  (V  =  15  m/s  and  PS  =  5%). 


■  MBDyn 

■  MBDyn 

•  mapped  DIC 

•  mapped  DIC 

(a)  AoA  =  6  deg.  (b)  AoA  =  9  deg. 


Figure  2.15:  DIC/MBDyn  correlation  of  wind  tunnel  tests  (V  =  18  m/s  and  PS  =  5%). 
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2.3.  Wind  Tunnel  Tests 


2.3.3  Inverse  Survey 

The  deformed  configuration  of  the  previously  investigated  rectangular  membrane  is  determined 
herein  using  the  strain  measurements  derived  from  DIC,  re-sampled  onto  the  numerical  mesh 
using  the  previously  discussed  MLS  procedure.  The  membrane  dimensions,  loading,  boundary 
conditions,  material  properties  are  the  same  as  in  the  direct  analysis  case  discussed  above,  see 
Table  2.1.  A  triangular  mesh  consisting  of  8x16  elements  has  been  used  for  the  IFEM  procedure 
(Figure  2.6):  the  strain  measurements  from  DIC  are  re-mapped  (and  smoothed  out)  onto  this 
virtual  strain  measurement  grid,  and  used  as  inputs  for  the  IFEM  analysis.  Figures  from  2.16 
to  2.19  compare  the  deformation  shape  corresponding  to  the  inverse  FEM  analysis  with  the 
experimental  data  re-sampled  on  the  numerical  mesh. 


•  mapped  DIC 

•  mapped  DIC 

♦  iFEM 

♦  iFEM 

(a)  AoA  =  3  deg. 


(b)  AoA  =  6  deg. 


Figure  2.16:  IFEM/DIC  correlation  of  wind  tunnel  tests  (V  =  12  m/s  and  PS  =  2%). 


•  mapped  DIC 

•  mapped  DIC 

♦  iFEM 

♦  iFEM 

(a)  AoA  =  3  deg.  (b)  AoA  =  6  deg. 

Figure  2.17:  IFEM/DIC  correlation  of  wind  tunnel  tests  (V  =  15  m/s  and  PS  =  2%). 
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•  mapped  DIC 

•  mapped  DIC 

♦  iFEM 

♦  iFEM 

(a)  AoA  =  6  deg. 

Figure  2.18:  IFEM/DIC  correlation  of  wind  tunnel  tests  (V 


(b)  AoA  =  9  deg. 


=  15  m/s  and  PS 


5%). 


•  mapped  DIC 

•  mapped  DIC 

♦  iFEM 

♦  iFEM 

(a)  AoA  =  6  deg.  (b)  AoA  =  9  deg. 

Figure  2.19:  IFEM/DIC  correlation  of  wind  tunnel  tests  (V  =  18  m/s  and  PS  =  5%). 
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2.3.  Wind  Tunnel  Tests 


Utilizing  DIC  data,  the  prediction  of  max  camber,  where  c  is  the  chord,  can  be  evaluated 
for  its  accuracy  and  overall  physical  behavior.  Figures  2.20  and  2.21  represent  the  average  max 
(mapped)  measured  ( “EXP” )  and  predicted  ( “FSI”  and  “IFEM” )  static  camber  for  two  different 
values  of  prestrain,  2%  and  5%,  respectively,  and  for  a  flow  velocity  of  V  =  12,  15  and  18  m/s. 
The  measured  and  the  predicted  displacement  are,  as  expected,  primarily  characterized  by  an 
adaptive  inflation,  which  increases  the  local  camber. 


(a)  V  =  12  m/s.  (b)  V  =  15  m/s. 

Figure  2.20:  Measured  max  displacement  from  DIC  data  and  predicted  max  displacement  (PS  =  2%). 


(a)  V  =  15  m/s.  (b)  V  =  18  m/s. 

Figure  2.21:  Measured  max  displacement  from  DIC  data  and  predicted  max  displacement  (PS  =  5%). 
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2.3.  Wind  Tunnel  Tests 


Figures  2.22  and  2.23  show  the  coefficient  of  lift  measured  by  the  load  cell  attached  to  the 
wing  (“EXP”).  They  show  the  coefficient  of  lift  calculated  by  the  coupled  fluid-structure  simula¬ 
tion  ( “FSI” ) ,  and  also  the  estimated  coefficients  of  lift  found  by  integrating  the  estimated  pressure 
distributions  (calculated  from  the  remapped  full-field  DIC  measurements,  “mapDIC,”  and  from 
the  estimated  displacements  via  inverse  analysis,  “IFEM,”  respectively)  to  find  the  normal  load, 
converting  it  into  lift  via  the  AoA,  and  finally  calculating  the  lift  coefficient. 

It  should  be  noted  that  those  curves  are  nearly  invariant  with  respect  to  changes  in  flow 
velocity  within  the  range  of  this  study  (Re  =  45k-67k).  The  error  bars  in  the  previous  figures 
represent  the  standard  deviations,  crcexp,  which  are  presented  for  reference  in  Table  2.4. 


(a)  V  =  12  m/s.  (b)  V  =  15  m/s. 

Figure  2.22:  Static  lift  model  and  wind  tunnel  data  (PS  =  2%). 


(a)  V  =  15  m/s.  (b)  V  =  18  m/s. 

Figure  2.23:  Static  lift  model  and  wind  tunnel  data  (PS  =  5%). 
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Table  2.4:  Wind  tunnel  test,  results. 
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Chapter  2.  Results  2.3.  Wind  Tunnel  Tests 


Conclusions 


This  work  presents  the  direct  and  inverse  analysis  of  membrane  elements  for  fluid-structure  in¬ 
teraction  problems. 

To  summarize:  (i)  a  membrane  inverse  analysis  based  on  a  three-node  membrane  element 
was  developed,  based  on  a  least  squares  smoothing  functional  that  employs  the  complete  set  of 
strain  measures;  (ii)  a  four-node  membrane  element  was  implemented  in  a  multibody-based  co¬ 
simulation  analysis  for  the  direct  simulation  of  coupled  fluid-structure  problems;  (iii)  the  inverse 
analysis  has  been  verified  by  reconstructing  the  deformed  solution  obtained  with  the  analogous 
direct  formulation  applied  on  a  different  mesh  and  subsequently  re-sampled;  (iv)  both  the  direct 
and  the  inverse  analyses  have  been  validated  by  comparing  the  direct  prediction  and  the  recon¬ 
structed  deformation  with  experimental  data  for  prestressed  rectangular  membranes  subjected  to 
hydrostatic  pressure  loads;  (v)  an  approach  to  estimating  aerodynamic  load  present  on  a  flexible 
membrane  wing  from  elastic  strain  sensor  was  developed. 

The  proposed  analysis  enables  accurate  and  computationally  efficient  high-fidelity  deformation 
reconstruction  solutions.  It  is  therefore  applicable  to  both  static  and  dynamic  problems.  Hydro¬ 
static  pressure  tests  were  considered  in  order  to  compare  a  known  input  load  to  the  estimated 
load.  Results  were  favorable:  the  average  hydrostatic  pressure  estimate  was  reasonably  close  to 
the  actual  applied  hydrostatic  pressure.  In  addition,  the  error  introduced  in  the  estimated  pres¬ 
sure  distribution  can  be  seen  in  the  irregularity  of  the  estimated  pressure  distribution  compared 
to  the  ideal  pressure  distribution.  The  proposed  procedure  for  the  reconstruction  of  shape  and 
distributed  loads  is  able  to  operate  at  sample  rates  of  the  order  of  30  Hz,  thus  meeting  the  initial 
real-time  operation  requirement. 
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Appendix  A 


Inverse  Analysis:  2D  Example 


Consider  a  two-dimensional  problem  in  the  x-z  plane.  The  strain  is 


U2/ 

/x 


l2,  - 

/* 


W 


2 

/x 


2 


where  the  quadratic  term  in  u/x  is  neglected  owing  to  the  presence  of  the  corresponding  linear 
term.  The  error  is  e  =  ex  -  .  The  cost  function  becomes: 


[e  H-  ^  (w/x  ^ref/x)  j 

f  (m)\2  i  f  (m)\  2 

[u/x  -  eyx  'J  +\u/x-e]D  )w/x 


~^+k  (w/x  -  i 


ef  /  x ) 


(A.l) 


Its  gradient  is: 
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(A. 2a) 


(A. 2b) 


Clearly  u/x  =  e£m-)  and  w/x  =  0  makes  the  gradient  vanish  when  k  =  0,  although  it  likely  is  not  the 
solution  that  is  sought. 

Consider  now  the  Hessian  matrix;  its  elements  are 


(f>  ,  =1 

^/u/xu/x 
^/u/xw/x  ^/x 

i  (m)  .  3  o  .j 

^ /w JXW Jx  U/x  £x  +  2^/x  T  ^ 

For  the  Hessian  matrix  to  be  positive  definite,  /c  must  be: 


(A. 3a) 
(A. 3b) 
(A. 3c) 


Consider  now  a  (simplified  and  heuristic)  iterative  solution  procedure: 

1.  initialize  u/x  =  0,  w/x  =  0; 

2.  initialize  iuref/.,.  ^  0  (“similar”  to  the  expected  solution); 
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3.  after  linearization,  the  increment  of  the  derivative  of  the  transverse  displacement  is 


Aw/x  = 


,  X  W  , 

k  +  u/x-  e|bm)  + 


(^ref/x  ^/x) 


k  +  e 


(  V’ref /.T  '>l-' / x  ) 


(A.4) 


which  indicates  that  w/x  progresses  towards  the  reference  solution  with  a  pace  that  is  unit 
when  k  ->■  +oo,  namely  Aw,x  =  wlet/x—w/x',  for  k  0,  the  increase  tends  to  vanish  for  positive 

and,  for  0  <  k  <  -  (u/x  -  =  -e,  could  even  be  in  the  wrong  direction  (in 


error 


such  case,  the  Hessian  matrix  would  be  indefinite); 
4.  the  increment  of  the  in-plane  displacement  is 


Au/x  =  ~ 


(m)  i  /a 

U/X~£x  +  — 


(  \  i  (m)  i  !x 

W / x  {^rei/x  ^  /  x)  u/x  ^ x  ^ 


(m)  .  /x 
U/x-Zx  +  — 


k  [ w/x  ( wief/x  -  w/x)  +e]+e2 


k  +  e 


For  k  — >  +00  the  increment  is 


(A.5) 


(m)  W/x 

A u/x  [  U}  /  X(/Wre{/X  w/x)  d  U/x  ^x  d  “  1  (w/xi^ref/x  ^ / x)  d  e)  5 


for  fc  ^  0  it  is 


A  I  (m)  i  W/x 

Au/x  =  -  \u/x  -  ei  d  — 


The  previous  consideration  applies  about  the  sign  of  the  increment. 


/ u2  1^2  \ 

Consider  the  strain  ex  =  u/x  +  (  “/x  2W/x  j  and  equate  it  to  its  measured  value,  ex  =  e£m). 

Assume  that  w/x  =  wlei/x,  regardless  of  its  value.  One  obtains: 

U%  +  2U/X  ~  (2e^m)  -  Wret/x)  =  0  (A-6) 

whose  solution  is 


and 


u/x  =  -1  ±  \/l  +  2e(xm)  -  w2e{/x 


(A.7) 


l(x)  =  Io  (_1±/ 


ld24m)(0-<f/a:(0  )  d£ 


(A.8) 


Nothing  prevents  the  solution  from  jumping  between  the  ±  cases  as  needed  to  account  for  the 
boundary  condition  u{t)  =  0.  Note  that  two  extreme  cases  could  be  wlef/x  =  0,  and  wlei/x  =  ±\J 2e).m\ 
The  latter  yields  u(x)  =  0. 


In  conclusion,  a  single-curvature  solution  cannot  be  reconstructed. 
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